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maj elloOslac . Stanford . edu 
ABSTRACT 

This is the first of a series of papers aimed at characterizing the populations 
detected in the high-latitude sky of the Fermi-LAT survey. In this work we 
focus on the intrinsic spectral and flux properties of the source sample. We 
show that when selection effects are properly taken into account, Fermi sources 
are on average steeper than previously found (e.g. in the bright source list) 
with an average photon index of 2.40±0.02 over the entire 0.1-100 GeV energy 
band. We confirm that FSRQs have steeper spectra than BL Lac objects with 
an average index of 2.48±0.02 versus 2.18±0.02. Using several methods we 
build the deepest source count distribution at GeV energies deriving that the 
intrinsic source (i.e. blazar) surface density at Fi o > 10~ 9 ph cm -2 s -1 is 
0.12^02 deg~ 2 . The integration of the source count distribution yields that 
point sources contribute 16(±1.8) % (±7% systematic uncertainty) of the GeV 
isotropic diffuse background. At the fluxes currently reached by LAT we can 
rule out the hypothesis that point-like sources (i.e. blazars) produce a larger 
fraction of the diffuse emission. 

Subject headings: cosmology: observations - diffuse radiation - galaxies: active 
gamma rays: diffuse background - surveys - galaxies: jets 

1. Introduction 

The origin of the extragalactic gamma-ray background (EGB) at GeV 7-rays is one of 
the fundamental unsolved problems in astrophysics. The EGB was first detected by the SAS- 
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2 mission (IFichtel et al.lll975l ) and its spectrum was measured with good accuracy by the 



Ener getic Gamma Ray Experiment Telescope (EGRET ISreekumar et al.lll998l ; IStrong et al. 



2004 ) on board the Compton Observatory. These observations by themselves do not provide 
much insight into the sources of the EGB. 

Blazars, active galactic nuclei (AGN) with a relativistic jet point ing close to our line o f 
sight, represent the most numerous population detected by EGRET lHartman et al.l (Il999l ) 
and their flux constitutes 15 % of the total EGB intensity (resolved sources plus diffuse 
emission). Therefore, undetected blazars (e.g. all the blazars under the sensitivity level of 
EGRET) are the most likely candidates for the origin of the bulk of the EGB emission. 
Studies of the luminosity function of blazars showed that the c ontribution of blazars to 
the EGRET EGB could be in the range from 20 % to 100 % (e.g. Istecker fc Salamonlll996 



Chiang fc Mukherjed Il998l ; iMiicke &: Pohll |2000| ) . although the newest derivations suggest 



that blazars are responsible for only 
Dermerll20~0~7i llnoue fc Totanill2009h . 



20-40 % of the EGB (e.g. iNarumoto fc Totanil 12006 



It is thus possible that the EGB emission encrypts in itself the signature of some of 



duced by the assembly of Large Scale S tructures (e.g. lLoeb fc Waxmanl 12000c iMiniati 
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Thompson et al.ll2007f ). are among the most likely candidates for the generation of dif- 
fuse the GeV emission. Dark matter (DM) which constitutes more than 80 % of the matter 
in the Universe can also provide a diffuse, cosmological, background of 7-rays. Indeed, super- 
symmetric theories with R-parity predict that the lig htest DM particles (i.e . , the neutralinos) 



Ullio et al. 


2002; 


Ahn et al. 


2007) 



With the advent of the Fermi Large Area Telescope (LAT) a better understanding of the 
origin of the GeV diffuse emission becomes possible. Fermi has recently perfo rmed a new 



measu rement of the EGB spectrum (also called isotropic diffuse background, lAbdo et al. 
2010dl ). This has been found to be consistent with a featureless power law with a pho- 
ton index of ~2.4 in the 0.2-100 GeV energy range. The integrated flux (E>100MeV) of 
1.03(±0.17) x 10~ 5 ph cm~ 2 s _1 sr -1 has been found to be significantly lo wer than the one 
of 1.4 5(±0.05)xlQ- 5 ph cm 
1998h . 



2 s 1 sr 1 determined from EGRET data 



sec 



Sreekumar et al. 



In this study we address the contribution of unresolved point sources to the GeV dif- 
fuse emission and we discuss the implication s. Early findings on the integrated emission 
of unresolved blazars were already reported in lAbdo et al.l (j2009al ) using a sample of bright 
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AGN detected in the first three months of Fermi observations. The present work represents 
a large advance, with ~4 times more blazars and a detailed investigation of selection effects 
in source detection. 

This work is organized as follows. In § [3] the intrinsic spectral properties of the Fermi 
sources are determined. In § [4] the Monte Carlo simulations used for this analyses are 
outlined with the inherent systematic uncertainties (see § [5j). Finally the source counts 
distributions are derived in § and § [7| while the contribution of point sources to the GeV 
diffuse background is determined in § [SI § M discusses and summarizes our findings. Since 
the final goal of this work is deriving the contribution of sources to the EGB, we will only use 
physical quantities (i.e. source flux and photon index) averaged over the time (11 m onths) 



included in the analysis for the First Fermi-LAT catalog (lFGL. lAbdo et al.ll2010bl ). 



2. Terminology 

Throughout this paper we use a few terms which might not be familiar to the reader. 
In this section meanings of the most often used are clarified. 

• spectral bias: (or photon index bias) is the selection effect which allows Fermi-LAT to 
detect spectrally hard sources at fluxes generally fainter than for soft sources. 

• flux-limited sample: it refers to a sample which is selected uniformly according solely 
to the source flux. If the flux limit is chosen to be bright enough (as in the case of 
this paper), then the selection effects affecting any other properties (e.g. the source 
spectrum) of the sample are negligible. This is a truly uniformly selected sample. 

• diffuse emission from unresolved point sources: it represents a measurement of the 
integrated emission from sources which have not been detected by Fermi. As it will 
be shown in the next sections, for each source detected at low fluxes, there is a large 
number of sources which have not been detected because of selection effects (e.g. the 
local background was too large or the photon index was too soft, or a combination 
of both). The diffuse emission from unresolved point sources (computed in this work) 
addresses the contribution of all those sources which have not been detected because 
of these selection effects, but have a flux which is formally larger than the faintest 
detected source. 
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3. Average Spectral Properties 



3.1. Intrinsic Photon index distributions 



As shown already in (lAbdo et al.ll2009al . but see also Fig. [T]), at faint fluxes the LAT 
detects more easily hard-spectrum sources rather than sources with a soft spectrum. Sources 
with a photon index (e.g. the exponent of the power-law fit to the source photon spectrum) 
of 1.5 can be detected to fluxes which are a factor > 20 fainter than those at which a source 



with a photon index of 3.0 can be detected (see lAbdo et al.l l2010d . for details). Thus, given 
this strong selection effect, the intrinsic photon index distribution is necessarily different 
from the observed one. An approach to recover the intrinsic photon index distribution is 
obtained by studying the sample above Fi o ~ 7 x 10~ 8 ph cm" 2 s" 1 and |6| > 10° (see right 
panel of Fig. [T]). Indeed above this flux limit, LAT detects all sources irrespective of their 
photon index, flux or position in the high-latitude sky. Above this limit LAT detects 135 
sources. Their photon index distribution, reported in Fig. [1] is compatible with a Gaussian 
distribution with mean of 2.40±0.02 and dispersion of 0.24±0.02. These values differ from 
the mean of 2.23±0.01 and dispersion of 0.33±0.01 derived using the entire |6| > 10° sample. 
Similarly the intrinsic photon-index distributions of FSRQs and BL Lacs are different from 
the observed distributions. In both case the observed average photon-index is harder than 
the intrinsic average value. The results are summarized in Tab. [TJ 




Fig. 1. — Left PaneLFlux-photon index plane for all the |6| > 10° sources wi th TS> 25. 



The da shed line is the flux limit as a function of photon index reported in lAbdo et al. 



( l2010el ). while the solid line represents the limiting flux above which the spectral selection 
effects become negligible. Right Panel: Photon index distribution of all sources for Fioo > 
7 x 10~ 8 ph cm~ 2 s _1 . Above this limit the LAT selection effect towards hard sources becomes 
negligible. 
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3.2. Stacking Analysis 

Another way to determine the average s pectral properties i s by stacking source spectra 



together. This is particularly simple since ( lAbdo et al.l l2010b[ ) reports the source flux in 



five different energy bands. We thus performed a stacking analysis of those sources with 
Fioo > 7x 10~ 8 ph cm" 2 s _1 , TS> 25, and |6| >10°. For each energy band the average flux 
is computed as the weighted average of all source fluxes in that band using the inverse of 
the flux variance as a weight. The average spectrum is shown in Fig. [2j A power law model 
gives a satisfactory fit to the data (e.g. \ 2 /dof ~ 1), yielding a photon index of 2.41±0.07 
in agreement with the results of the previous section. 

We repeated the same exercise separately for sources identified as FSRQs and BL Lacs 
in the flux-limited sample. Both classes have an average spectrum which is compatible with 
a single power law over the whole energy band. FSRQs are characterized by an index of 
2.45±0.03 while BL Lac objects have an average index of 2.23±0.03 



Monte Carlo Simulations 



In order to estimate the LAT sky coverage robustly we performed detailed Monte Carlo 
simulations. The sche me of the simulation procedure is an improved version of what has 
already been applied in lAbdo et al.l ( I2009af ). We performed 18 end-to-end simulations of the 
LAT sky which resemble as closely as possible the observed one. The tool gfo6ssm£]has been 
used for this purpose. For each simulation we modeled the Galactic and isotropic diffuse 



1 The list of science tools for the analysis of Fermi data is accessible at 
http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.litml 

Table 1. Observed versus Intrinsic photon indices distributions for the Fermi/LAT source 

classes. 



Observed Distribution Intrinsic Distribution 



SAMPLE mean a mean a 



ALL 2.24±0.01 0.31±0.01 2.40±0.02 0.24±0.02 

FSRQ 2.44±0.01 0.21±0.01 2.47±0.02 0.19±0.02 
BL Lac 2.05±0.02 0.29±0.01 2.20±0.04 0.22±0.03 
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Fig. 2. — Stacked spectrum of sources in the flux-limited sample. The dashed line is the best 
power law fit with a slope of 2.41±0.07. 




Fig. 3. — Stacked spectrum of FSRQs (left) and BL Lac objects (right) in the Fermi-LAT 
flux-limited sample. 

backgrounds using models (e.g. gll_iem_v02.fit) currently recommended by the LAT team. 

An isotropic population of point-like sources was added to each simulated observation. 
The coordinates of each source were randomly drawn in order to produce an isotropic dis- 
tribution on the sky. Source fluxes were randomly drawn from a standard log iV-log S 
distribution with parameters similar to the one observed by LAT (see next section). Even 
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though the method we adopt to derive the survey sensitivity does not depend on the nor- 
malization or the slope of the input log AMog S, using the real distribution allows simulated 
observations to be produced that closely resemble the sky as observed with the LAT. The 
photon index of each source was also drawn from a Gaussian distribution with mean of 2.40 
and lcr width of 0.28. As noted in the previous section, this distribution represents well 
the intrinsic (not the observed one) distribution of photon indices. The adopted dispersion 
is slightly larger than what was found in the previous section and it is derived from the 
analysis of the entire sample (see § 16. 2p . In this framework we are neglecting any possible 
dependence of the photon index distribution with flux. Also we remark that the approach 
used here to derive the source count distribution depends very weakly on the assumptions 
(e.g. the log AMog S used) made in the simulations. 

More than 45000 randomly distributed sources have been generated for eac h realization 



of the s imulations. Detection follows (albeit in a simpler way) the scheme used in lAbdo et al. 



( l2010bl ). This scheme adopts three energy bands for source detection. The first band includes 



all /ront-converting^] and foac^-converting photons with energies larger than 200 MeV and 
400 MeV, respectively. The second band starts at 1 GeV for front photons and at 2 GeV 
for back photons. The high-energy band starts at 5 GeV for front photons and at 10 GeV 
for back photons. The choice of combining front and back events with different energies is 
motivated by the fact that front events have a better point spread function (PSF) than back 
ones. The two PSFs are similar when the energy of ftacfc-converting photons is approximately 
twice as that of /roni-converting ones. The image pixel sizes changes according to the energy 
band and is 0.1, 0.2 and 0.3 degrees for the low, medium and high-energy bands respectively. 
The final list of candidate sources is obtained starting the detection at the highest energy 
band and adding all those sources which, being detected at lower energy, have a position not 
consistent with those detected at high energy. The detection step uses pgwave fo r determining 



the p osition of the excesses and pointfit for refining the source position. Pgwave (ICiprini et al. 



20071 ) is a tool which uses several approaches (e.g. wavelets, threshold ing, image denoising 



and a sliding cell algorithm) to find source candidates while pointfit (e.g. iBurnett et al.ll2009f ) 



employes a simplified binned likelihood algorithm to optimize the source position. 

All the source candidates found at this stage are then ingested to the Maximum Like- 
lihood (ML) algorithm gtlike to determine the significance and the spectral parameters. In 
this step all sources' spectra are modeled as single power laws. On average, for each simu- 
lation only ~1000 sources are detected (out of the 45000 simulated ones) above a TS| of 25 



2 Photons pair-converting in the top 12 layers of the tracker are classified as fro rat-converting photons or 
fracfc-converting otherwise. 

3 The test statistics (or TS) is defined as: TS=— 2(lnL — InLi). Where L and Li are the likelihoods of 



-11- 



and this is found to be in good agreement with the real data. 



4.1. Performances of the detection algorithm on real data 

In order to test the reliability of our detection pipeline we applied it t he to real 1 



year da taset. Our aim was to cross check our result with the result reported in lAbdo et al. 



f l2010bf ). The flux abo ve 100 MeV, comput ed from the power-law fit to the 100 MeV-100 GeV 
data, is not reported in lAbdo et al.l ( l2010bl ) . but it can be readily obtained using the following 
expression: 

100 



100 



E 



piv 



l-r 



x i-r 



-i 



i 



where F 100 is the 100 MeV-100 GeV photon flux, T is the absolute value of the p hoton index, 
E piv is the pivot energy and F density is the flux density at the pivot energy (see lAbdo et al. 



201 Obi , for details). Fig. H] shows the comparisons of both fluxes (above 100 MeV) and of 
the photon indices for the sources detected in both pipelines. It is clear that the fluxes and 
photon i ndices deriv e d in th is analysis are reliable; for each source they are consistent with 
those in lAbdo et al.l ( 12 01 Obi ) within the reported errors. 



The number o f sources detected in the simplified pipeline is smaller than found by 
Abdo et al.l ( 1201 Obi ). Above a TS of 50 and \b\ > 20° our approach detects 425 sources while 
the 1FGL catalog has 497. Indeed, our aim is n o t to pr oduce a detection algorithm which 
is as sensitive than the one used in I Abdo et al.l ( )2010bl ). but a detection algorithm which 
is proven to be reliable and can be applied consistently to both real data and simulations. 
This allows us to assess properly all selection effects important for the LAT survey and its 
analysis. O n this note we rema rk that all the 425 sources detected by our pipeline are also 
detected by lAbdo et al.l (l2010bl ) . For this reason we limit the studies presented in this work 
to the subsample of sources which is detected by our pipeline. The det ails of this sample o f 
sour ces are reported in T ab. [2j The associations are the ones reported in lAbdo et al.l (|2010eJ) 
and I Abdo et al.l (]2010bl ). In our sample 161 sources are classified as FSRQs and 163 as BL 
Lac objects while only 4 as blazars of uncertain classification. The number of sources which 
are unassociated is 56, thus the identification incompleteness of this sample is ~13%. 



th e background (null hypothesis) and the hypothesis being tested (e.g. source plus bac kground). Accordi ng 
to Wilksl ( 19381 ) . the significance of a detection is approximately n a = ^J\TS) (see also Aiello et al. 20081 ). 
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Table 2. Composition of the |b| >20 TS>50 sample used in this analysis. 



CLASS 


# objects 


Total 


425 


FSRQs 


161 


BL Lacs 


163 


Uncertain 51 


4 


Blazar Candidates 


24 


Radio Galaxies 


2 


Pulsars 


9 


Others' 3 


6 


Unassociated sources 


56 



a Blazars with uncertain classifica- 
tion. 

b It includes Starburst galaxies, 
Narrow line Seyfert 1 objects and 
Seyfert galaxy candidates. 
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10= 10' 7 10 6 1 1 - 5 2 25 

F 100 (LAT Catalog) [ph cm' 2 s 1 ] Photon Index (LAT Catalog) 



Fig. 4. — Performance of the detect i on pip eline used in this work with respect to the de- 



tection pipeline used in lAbdo et al.l (l2010bT ) . The left panel shows the comparison of the 
reconstructed 7-ray fluxes while the right panel shows the comparison of the photon indices. 
In both cases the solid line shows the the locus of points for which the quantity reported in 
the y-axis equals the one in the x-axis. 

4.2. Derivation of the Sky Coverage 

In order to derive the sky coverage from simulations, detected sources (output) need 
to be associated to the simulated ones (input). We do this on a statistical basis using an 
estimator which is defined for each set of input-output sources as: 

where x, S and T are the source coordinates, fluxes and photon indices as determined from 
the ML step while x , S and T are the simulated (input) values. The 1 a errors on the 
position, flux and photon index are a pos , as and a? respectively. We then flagged as the 
most likely associations those pairs with the minimum value of R 2 . All pairs with an angular 
separation which is larger than the 4 a error radius are nagged as spurious and excised from 
the following analysis. The empirical, as der ived from the real da ta, 5 a error radius as a 



functi on of source TS is shown in Fig. |5j As in lHasinger et al.l (119931 ) and in lCappelluti et al. 



( 120071 ) we defined confused sources for which the ratio S/ (S + 3a s) (where as is the error 
on the output flux) is larger than 1.5. We found that, according to this criterion, ~4% of 
the sources (detected for |6| > 10°) are confused in the first year survey 

The right panel of Fig. [6] shows the ratio of reconstructed to simulated source flux versus 
the simulated source flux. At medium to bright fluxes the distribution of the ratio is centered 
on unity showing that there are no systematic errors in the flux measurement. At low fluxes 
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(in particular for Fioo < 10 9 ph cm 2 s l ) the distribution is slightly (or somewhat) biased 
toward values great er than unity. T his is produced by three effects: 1) source confusion, 



2) Eddington bias ( lEddingtonl Il940l ) and 3) non converging Maximum Likelihood fits (see 
§ I5.1l for details). The Eddington bias arises from measurement errors of any intrinsic source 
property (e.g. source flux). Given its nature, it affects only sources close to the detection 
threshold. Indeed, at the detection threshold the uncertainty in the reconstructed fluxes 
makes sources with a measured flux slightly larger than the real value more easily detectable 
in the survey rather than those with a measured flux slightly lower than the real one. This 
causes the shift of the flux ratio distribution of Fig. [6] to move systematically to values larger 
than unity at low fluxes. In any case, the effect of this bias is not relevant as it affects less 
than 1 % of the entire population. This uncertainty will be neglected as only sources with 
Fioo — 10~ 9 ph cm -2 s -1 will be considered for the analysis presented here. Moreover, the 
right panel of Fig. [H] shows that the measured photon index agrees well with the simulated 
one. 

In addition to assessing the reliability and biases of our source detection procedure, the 
main aim of these simulations is to provide a precise estimate of the completeness function 
of the Fermi/LAT survey (known also as sky coverage). The one-dimensional sky coverage 
can be derived for each bin of flux as the ratio between the number of detected sources 
and the number of simulated sources. The detection efficiency for the entire TS> 50 and 
\b\ > 20° sample is reported in Fig. [71 This plot shows that the LAT sensitivity extends all 
the way to F 10 o ~ 10 -10 ph cm -2 s" 1 although at those fluxes only the hardest sources can be 
detected. Also the sample becomes complete for F 10 o = 7 — 8 x 10 _8 ph cm -2 s _1 . Since for 
these simulations, the intrinsic distribution of photon indices has been used (see e.g. § 13.11) 
this sky coverage takes properly into account the bias towards the detection of hard sources. 
This also means that this sky coverage cannot be applied to other source samples with very 
different photon index distributions. 



5. Systematic Uncertainties 

5.1. Non converging Maximum Likelihood fits 

A small number of sources detected by our pipeline have unreliable spectral fits. Most 
of the time, these sources have a reconstructed photon index which is very soft (e.g. ~5.0) 
and at the limit of the accepted range of values. As a consequence their reconstructed flux 
overestimates the true flux by up to factor 1000 (see left panel of Fig. E]). This is due to 
the fact the the ML algorithm does not find an absolute minimum of the fitting function for 
these cases. Inspection of the regions of interests (ROIs) of these objects shows that this 
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Fig. 5. — Angular separation of the real LAT sources from the most probable associated 
counterpart as a function of TS . All sources with |6| >10° with a probability of associations 
larger than 0.5 were used (see lAbdo et al.ll2010eL for a definition of probability of associa- 
tion). The solid line is the best fit for the mean offset of the angular separations while the 
dashed line represents the observed 5 a error radius as a function of test statistics. Note that 
the 5 a error radius is weakly dependent on the level of probability of association chosen. 



tends to happen either in regions very dense with sources or close to the Galactic plane, 
where the diffuse emission is the brightest. The best approach in this case would be to adopt 
an iterative procedure for deriving the best-fitting parameters which starts by optimizing the 
most intense components (e.g. diffuse emissions and b right sources ) and t hen move to the 
fainter ones. This procedure is correctly implemented in lAbdo et al.l (l2010bl ). Its application 
to our problem would make the processing time of our simulations very long and we note 
that the systematic uncertainty deriving from it is small. Indeed, the number of sources with 
unreliable spectral parameters are for TS > 25 are 2.3% and 2.0% for \b\ > 15° a |6| > 20° 
respectively. These fractions decrease to 1.2% and 0.9% adopting TS> 50. 

To limit the systematic uncertainties in this analysis, we will thus select only those 
sources which are detected above TS> 50 and |6| > 20°. It will also be shown that results 
do not change if the sample is enlarged to include all sources with |6| > 15°. 
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Fig. 6. — Left Panel: Reconstructed versus Simulated fluxes for all sources with TS>50 and 
|6| >20°. For the analysis reported here only sources with F 10 o > 10~ 9 ph cm -2 s^ 1 are 
considered. Right Panel: Reconstructed versus Simulated photon indices for all sources with 
TS>50 and \b\ >20°. 




Fig. 7. — Detection efficiency as a function of measured source flux for |6| > 20°, TS> 50 
and for a sample of sources with a mean photon index of 2.40 and dispersion of 0.28. The 
error bars represent statistical uncertainties from the counting statistic of our Monte Carlo 
simulations. 
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5.2. Variability 



It is well known that blazars are inherently variable objects with variability in flux of 
up to a factor 10 or more. Throughout this work only average quantities (i.e. mean flux 
and mean photon index) are used. This is correct in the context of the determination of 
the mean energy release in the Universe of each source. Adopting the peak flux (i.e. the 
brightest flux displayed by each single source) would produce the net effect of overestimating 
the true intrinsic source density at any flux (see the examples in lReimer fc Thompsonll200ll ) 
with the result of overestimating the contribution of sources to the diffuse background. 

It is not straightforward to determine how blazar variability affects the analysis pre- 
sented here. On timescales large enough (such as the one spanned by this analysis), the 
mean flux is a good estimator of the mean energy release of a source. This is not true 
anymore on short timescales (e.g. ~1 month) since the mean flux corresponds to the source 
flux at the moment of the observation. The continuous scanning of the 7-ray sky performed 
by Fermi allows to determine long-term variability with unprecedented accuracy. As shown 
already in Abdo et aljj 2009a ) the picture arising from Fermi is rather different from the one 
derived by EGRET ( iHartman et al.lll999l ). Indeed, the peak-to- mean flux ratio for Fermi 
sources is considerably smaller than for EGRET sources. For most of the Fermi sourc es this 
is just a factor 2, as is confirmed in the lyear sample (see Fig.10 in lAbdo et al.ll2010el ). This 
excludes the possibility that most of the sources are detected because of a single outburst 
which happened during the 1 1 months of observ ation and are undetected for the remaining 
time. Moreover, as shown in I Abdo et all (l2010d ) there is little or no variation of the photon 
index with flux. We thus believe that no large systematic uncertainties are derived from the 
use of average physical quantities and the total systematic uncertainty (see next section) will 
be slightly overestimated to accommodate possible uncertainties caused by variability. 



5.3. Non power law spectra 

It is well known that the spectra of blazars are complex and often show curvature when 
analyzed over a large waveband. In this case the approximation of their spectrum with a 
simple power law (in the 0.1-100 GeV band) might provide a poor estimate of their true flux. 
To estimate the uncertainties derived by this assumption we plotted for the extragalactic 
sample used here (e.g. TS> 50 and >20°) the source flux as derived from the power-law 
fit to the whole band versus the source flux as derived from the sum of the fluxes over the 5 



energy bands reported in I Abdo et al.l (]2010bt ) . This comparison is reported in Fig. [SJ From 



the figure it is apparent that the flux (F 10 o) derived from a power-law fit to the whole band 
overestimates slightly the true source flux. Analysis of the ratio between the power-law 
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flux and flux derived in 5 energy bands, shows that on average the Fioo flux overestimates 
the true source flux by ~8%. At very bright fluxes (e.g. F 10 o > 10~ 7 ph cm -2 s" 1 ) the 
overestimate reduces to ~5 %. For the analysis presented here we will thus assume that 
the total systematic uncertainty connected to the use of fluxes computed with a power-law 
fit over the broad 0.1-100 GeV band is 8 %. 

Considering also the uncertainties of the previous sections, we derive that the total 
systematic uncertainty for the sample used here (TS>50 and \b\ >20°) is ~10%. Since 
this uncertainty affects mostly the determination of the source flux it will be propagated by 
shifting in flux the sky coverage of FigJ7]by ±10 %. 




Fig. 8. — Source flux estimated with a power-law fit to the 0.1-10 GeV band versus the 
sum of the source fluxes derived in 5 contiguous energy bands (see lAbdo et all l2010bl for 
details). The solid line is the F banrfs =F 10 o relation. The spread at low fluxes arises from the 
difficulties of estimating the source flux in small energy bands. 



6. Source Counts Distributions 

The source counts distribution, commonly referred to as log N-\ogS or size distribution, 
is the cumulative number of sources N(> S) detected above a given flux S. In this section 
we apply several methods to derive the source count distribution of Fermi/LAT sources. We 
also remark that the catalog used for this analysis is the one described in § 14.11 (see also 
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Tab. E]). 



6.1. Standard Approach 

A standard way to derive the (differential) log iV-log 5* is through the following expres- 
sion: 

dS A S ^ a K ' 

i=l 1 

where N& s is the total number of detected sources with fluxes in the AS interval, and 
fli is the solid angle associated with the flux of the i t h source (i.e., the detection efficiency 
multiplied by the survey solid angle). We also note that formally iV is an areal density and 
should be expressed as dN/dVt. However for simplicity of notation the areal density will, 
throughout this paper, be expressed as N. For the |6| > 20° sample the geometric solid angle 
of the survey is 27143.6 deg 2 . In each flux bin, the final uncertainty is obtained by summing 
in quadrature the error on the number of sources and the systematic uncertainties described 
in § El 

Both the differential and the cumulative version of the source count distributions are 
reported in Fig. |HJ In order to parametrize the source count distribution we perform a \ 2 fit 
to the differential data using a broken power-law model of the type: 

^ = AS'^ S>S b 
db 

= AS~ Pl+h S-P 2 S <S b (4) 



where A is the normalization and S b is the flux break. The best-fit parameters are 
reported in Tab. El The log N- log S distribution of GeV sources shows a strong break 
{A/3 = f3 1 - 02 ~ 1-0) at Fioo = 6.97 (±0.13) x l(T 8 ph cm" 2 s -1 . At fluxes brighter than 
the break flux, the source count distribution is consistent with Euclidean (/3i = 2.5) while 
it is not at fainter fluxes. As Tab. |3] shows, these results do not change if the sample under 
study is enlarged to |6| >15°. 



6.2. A Global Fit 

Because of the spectral selection effect discussed in § 13.11 the sky coverage derived in 
14. 21 can be used only with samples which have a distribution of the photon indices similar to 
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Fig. 9. — Differential (left) and cumulative (right) log AMog S for all sources with TS> 50 
and \b\ > 20°. The dashed line is the best-fit broken power law model as reported in the 
text. 



the one used in the simulations (i.e. a Gaussian with mean and dispersion of 2.40 and 0.28). 
Here we aim at overcoming this limitation by implementing for the first time a novel, more 
formal, analysis to derive the source count distribution. We aim at describing the properties 
of the sample in terms of a distribution function of the following kind: 



where f(S) is the intrinsic flux distribution of sources and g(T) is the intrinsic distribu- 
tion of the photon indices. In this analysis, f{S) is modeled as a double power-law function 
as in Eq. |H The index distribution g(T) is modeled a Gaussian function: 

g (T) = e~^~ (6) 



where \i and a are respectively the mean and the dispersion of the Gaussian distribution. 
As it is clear from Eq. we made the hypothesis that the dN / dSdT function is factorizable in 
two separate distributions in flux and photon index. This is the most simple assumption that 
could be made and as it will be shown in the next sections it provides a good description of 
the data. Moreover, we emphasize, as already did in § HI that this analysis implicitly assumes 
that the photon index distribution does not change with flux. This will be discussed in more 
details in the next sections. 



This function is th en fitted to a l l data points using a Maximum Likelihood approach as 
described in Sec. 3.2 of lAjello et al.l (120091 ) . In this method, the Likelihood function can be 
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defined as: 



L = exp(-N exp ) J]A(5,,r,) (7) 

i=i 

with X(S, T) defined as: 

A(s,r) = J^n(s,r) (8) 

where fi(S', T) is the photon index dependent sky coverage and N ohs is the number of 
observed sources. This is generated from the same Monte Carlo simulation of § H] with the 
difference that this time the detection probability is computed for each bin of the photon- 
index-flux plane as the ratio between detected and simulated sources (in that bin). This 
produces a sky coverage which is function of both the source flux and photon index. 

The expected number of sources N exp can be computed as: 



N 

1 ' exp 



J dY J dSX(S,Y) (9) 



The maximum likelihood parameters are obtained by minimizing the function C(= 
-21nL): 

N ob3 

C = -2 HKSi, r*)) - 2N ln(N exp ) (10) 

i 

while their associated 1 a errors are computed by varying the parameter of interest, 
while the others are allowed to float, until an increment of AC=1 is achieved. T his gives an 



estimate of the 68 % confidence region for the parameter of interest (jAvnilll976l ). 



Once the dN/dSdT has been determined, the standard differential source count distri- 
bution can be readily derived as: 

dN _ f°° ^ dN 
dS 7„ 00 dSdT 



6.3. The Total Sample of Point Sources 



The results of the best-fit model for the entire sample of sources (for TS>50 and |6| >20°) 
are reported in Tab. HI Fig. [10] shows how well the best-fit model reproduces the observed 
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index and flux distributions. The x 2 test yields that the probabilities that the real distri- 
bution and the model line come from the same parent population are 0.98 and 0.97 for the 
photon index and flux distributions, respectively. In Fig. [IT] the source count distribution 
obtained here is compared to the one derived using the standard approach of § 16. 1( the good 
agreement is apparent. 

We also derived the source count distributions of all objects which are classified as 
blazars (or candidate blazars) in our sample. This includes 352 out of the 425 objects reported 
in Tab. |5J The number of sources that lack association is 56 and thus the incompleteness 
of the blazar sample is 56/425~ 13%. A reasonable and simple assumption is that the 
56 unassociated sources are distributed among the different source classes in a similar way 
as the associated portion of the sample (see Tab. [2]). This means that 46 out of the 56 
unassociated sources are likely to be blazars. As it is possible to notice both from the best- 
fit parameters of Tab. H] and from Fig. [12j there is very little difference between the source 
count distributions of the entire sample and the one of blazars. This confirms on a statistical 
basis that most of the 56 sources without association are likely to be blazars. It is also clear 
from Fig. [IUJ that the model (e.g. Eq. [5]) represents a satisfactory description of the data. 
This also implies that the intrinsic photon index distribution of blazars is compatible with 
a Gaussian distribution that does not change (at least dramatically) with source flux in the 
range of fluxes spanned by this analysis. A change in the average spectral properties of 
blazars with flux (and/or redshift) might be caused by the different cosmological evolutions 
of FSRQs and BL Lacs or by the spectral evolution of the two source classes with redshift. 
While it is something reasonable to expect, this effect is in the current dataset not observed. 
The luminosity function, which is left to a future paper, will allow us to investigate this 
effect in great detail. 



6.4. FSRQs 



For the classification of blazars as flat spectrum ra dio quasars ( FSRQs ) or BL Lacertae 
objects (BL Lacs) we use the same criteria adopted in lAbdo et al.l (l2009al). This classifica - 
tion relies on the c onven tional definition of BL Lac objects outlined in IStocke et al.l ( 1l99ll ). 
Urry fc Padovanil (119951 ). and iMarcha et al.l (119961 ) in which the equivalent width of the 
strongest optical emission line is <5A and the optical spectrum shows a Ca II H/K break 
ratio C<0.4. 



It is important to determine correctly the incompleteness of the sample when dealing 
with a sub-class of objects. Indeed, in the sample of Tabj2j 56 objects have no associations 
and 28 have either an uncertain or a tentative association with blazars. Thus the total 
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Fig. 10. — Distribution of photon indices (left) and fluxes (right) for the TS>50 and \b\ >20° 
sources. The dashed line is the best fit dN / dSdT model. Using the x 2 t est the probabilities 
that the data and the model line come from the same parent population are 0.98 and 0.97 
for the photon index and flux distribution respectively. 
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Fig. 11. — Comparison of log iV-log S of the whole sample of (TS>50 and |6| >20°) sources 
built with the standard method (green datapoints, see § 16. If) and the global fit method (red 
datapoints,see § 16 . 2 j) . 

incompleteness is 84/425 = ~19% when we refer to either FSRQs or BL Lac objects sepa- 
rately. The incompleteness levels of all the samples used here are for clarity reported also 
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Fig. 12. — Comparison between log iV-log S distributions of the whole sample of sources 
(solid circles) and blazars (open circles). The solid line are the respective best-fit models as 
reported in Tab. HI 



in Tab. HI Since we did not perform dedicated simulations for the FSRQ and the BL Lac 
classes, their source count distributions can be derived only with the method described in 

The best fit to the source counts (reported in Tab. H} is a double power-law model with 
a bright-end slope of 2.41±0.16 and faint-end slope 0.70±0.30. The log iV-log S relationship 
shows a break around Fioo =6.12(±1.30) x 10~ 8 ph cm -2 s _1 . The intrinsic distribution of 
the photon indices of FSRQs is found to be compatible with a Gaussian distribution with 
mean and dispersion of 2.48±0.02 and 0.18±0.01 in agreement with what found previously in 
Tab. [U The faint-end slope is noticeably flatter and this might be due to the fact that many 
of the unassociated sources below the break might be FSRQs. Fig. [TBI shows how the best-fit 
model reproduces the observed photon index and flux distributions. The \ 2 test indicates 
that the probability that the real distribution and the model line come from the same parent 
population is > 0.99 for both the photon index and flux distributions respectively. The 
left panel shows that the photon index distribution is not reproduced perfectly. This might 
be due to incompleteness or by the fact that the intrinsic distribution of photon indices is 
actually not Gaussian. However, a Kolmogorov-Smirnov (KS) test between the predicted 
and the observed distribution yields that both distributions have a probability of ~ 96 % of 
being drawn from the same parent population. Thus the current dataset is compatible with 
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the hypothesis that the intrinsic index distribution is Gaussian. The log AMog S of FSRQs 
is shown in Fig. [TH 




Fig. 13. — Distribution of photon indices (left) and fluxes (right) for the TS>50 and \b\ >20 
sources associated with FSRQs. 




Fig. 14.— Cumulative (left) and differential (right) source count distribution of Fermi 
blazars and the sub-samples reported in Tab. HJ Given the selection effect towards spectrally 
hard sources, BL Lac objects are detected to fluxes fainter than FSRQs. The flattening at low 
fluxes of the FSRQs logiV-logS is probably due to incompleteness (see text for details). The 
"All Blazars" class also includes all those sources which are classified as blazar candidates 
(see Tab. [2] for details). 
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6.5. BL Lacs 



The best-fit model of the source count distribution of the 161 BL Lac objects is again a 
broken power-law model. The break is found to be at Fioo =6.77±1.30 x 10 _8 ph cm -2 s _1 
while the slopes below and above the break are 1.72±0.14 and 2.74±0.30 respectively. The 
intrinsic photon index distribution is found to be compatible with a Gaussian distribution 
with mean and dispersion of 2.18±0.02 and 0.23±0.01 respectively. These results are in good 
agreement with the one reported in Tab. [TJ The best-fit parameters to the source counts 
distribution are reported in Tab. HI Fig. [15] shows how the best-fit model reproduces the 
observed photon index and flux distributions. The x 2 test indicates that the probability that 
the real distribution and the model line come from the same parent population is > 0.99 
for both the photon index and flux distributions respectively. The log iV-log S of BL Lacs, 
compared to the one of FSRQs and blazars, is shown in Fig. [Ill 




Fig. 15. — Distribution of photon indices (left) and fluxes (right) for the TS>50 and |6| >20° 
sources associated with BL Lacs. The dashed line is the best fit dN/dSdT model. 



6.6. Unassociated Sources 



We also constructed the log iV-log S of the 56 sources which are unassociated and it is 
reported in Fig. [HI Their source count distribution displays a very steep bright-end slope 
(/3i=3.16±0.50), a break around ~4.5xl0~ 8 ph cm -2 s _1 and faint-end slope of 1.63±0.24. 
The intrinsic photon index distribution is found to be compatible with a Gaussian distri- 
bution with mean and dispersion of 2.29±0.03 and 0.20±0.01 respectively (see Tab. H] for 
details). The extremely steep bright-end slope is caused by the fact that most (but not all) 
of the brightest sources have an association. Below the break the log N— log S behaves like 
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the one of blazars with the difference that the index distribution is suggesting that probably 
most of the sources are BL Lac objects. Indeed as can be seen in Fig. [TH all the sources with 
Fioo < 4 x 10 -8 ph cm -2 s _1 are identified as BL Lac objects in our sample. 



6.7. Unfolding Analysis 

Finally we employ a different approach to evaluate the log AMog S distribution based on 
a deconvolution (unfolding) technique. This method allows reconstructing the distribution 
of the number of sources from the data without assuming any model, also taking into account 
the finite resolution (i.e. dispersion) of the sky coverage. 

The purpose of the unfolding is to estimate the true distribution (cause) given the 
observed one (effect), and assuming some knowledge about the eventual migration effects 
(smearing matrix) as well as the efficiencies. The elements of the smearing matrix represent 
the probabilities to observe a given effect that falls in an observed bin Effectj from a 
cause in a given true bin Causei. In our case the observed distribution represents the 
number of sources as function of the observed flux above 100 MeV, while the true distribution 
represents the number of true sources as function of the true flux above 100 M eV. The 



unfolding algorithm adopted here is based on the Bayes theorem iD'Agostinil (119951 ) 



The smearing matrix is evaluated using the Monte Carlo simulation described in the 
§ 4. Its elements, P(F100j ]O b s |F100i itrne ), represent the probabilities that a source with a 
true flux above 100 MeV, F100i jtrue , is reconstructed with an observed flux above 100 MeV, 
FlOOj obg. The data are assumed to be binned in histograms. The bin widths and the number 
of bins can be chosen independently for the distribution of the observed and reconstructed 
variables. 

The log AMog S reconstructed with this method is shown in Fig. [16] and it is apparent 
that the source counts distributions derived with the 3 different methods are all in good 
agreement with each other. 



6.8. Comparison with Previous Estimates 

Fig. HH shows that the log AMog S distributions displays a strong break at fluxes 
Fioo ~ 6 x 10 -8 ph cm -2 s _1 . This represents the first time that such a flattening is seen 
in the log AMog S of 7-ray sources, blazar in particular. This is due to the fact that Fermi 
couples a good sensitivity to the all-sky coverage thus allowing to determine the source counts 
distribution over more than 3 decades in flux. 
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Above fluxes of Fi o = lCT 9 ph cm -2 s _1 , the surface density of sources is 0.12+q;q2 
deg -2 . At these f aint fl uxes our comparison can on ly be done with predictions from different 
models. iDermerl (120071 ) and llnoue &: Totanil (120091 ) predict a blazar surface density of respec- 
tively 0.030 deg -2 and 0.033 deg -2 . Both these predictions are a factor ~ 4 below the LAT 
measurement. However, it should be stressed that these models are based on the EGRET 
blazar sample which, because of strong selection effects against high-energy photons, counted 
a very limited number of BL Lac objects. 



ph cm 



-2 B -l> 



Dermerl (120071 ) predicts a density of 



At brighter fluxes (F 10 o > 5 x 10" 
FSRQs and BL Lacs of 4.1xl0~ 3 deg -2 and l.lxl0~ 3 deg~ 2 respectively. At the same flux, 
Mucke fc Pohll (120001 ) predict a density of 1.21 x 10~ 3 deg -2 and 3.04x 10~ 4 deg -2 respectively 
for FSRQs and BL Lac objects. The densities measured by Fermi are significantly larger, 
being 6.0(±0.6) x 10" 3 deg- 2 for FSRQs and 2.0(±0.3) x 10~ 3 deg- 2 for BL Lacs. 
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Fig. 16. — Source count distribution of Fermi point-like sources derived with three different 
methods. The distribution has been multiplied by (Fioo/10 -8 ) 1 ' 5 . The dashed line is the 
best fit model described in the text. The grey region indicates the flux at which a power law 
connecting the log A^-log S break (at ~ 6.6 x 10 -8 ph cm -2 s _1 ) and that given flux exceeds 
the EGB emission (see text for details). 
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7. Analysis in the Different Energy Bands 

The aim of the following analysis is to determine the contribution of point sources to the 
EGB in different contiguous energy bands. This is done by creating a logiV-logS distribution 
in 3 different energy bands: 0.1-1.0 GeV, 1.0-10.0 GeV and 10.0-100 GeV bands. This will 
allow us to study the spectrum of the unresolved emission from point sources and at the 
same time explore the properties of the source population in different bands. With this 
approach, the systematic uncertainty related to the flux estimate, given by the complex 
spectra of blazars (see $ 15.31) . will be removed. In addition, use of these bands should allow 
us to extend the survey region to \b\ > 10° (see § 15.11) . 

The analysis follows the method outlined in §H]with the difference that the final ML fit 
is restricted to the band under investigation. In the spectral fit, all parameters (including 
the photon index) are left free and are optimized by maximizing the likelihood function. 
Only sources that a given band have TS>25 are considered detected in that band. Formally 
each band and related sample is treated as independent here and no prior knowledge of the 
source spectral behaviour is assumed. In the three bands, the samples comprise respectively 
362, 597 and 200 sources detected for |b| >10° and TS> 25. 

In both the soft and the medium band (i.e. 0.1-1.0 GeV and 1.0-10.0 GeV), the logiV- 
\ogS is well described by a double power- law model, while in the hardest band (10-100 GeV) 
the logiV-logS is compatible with a single power-law model with a differential slope of 
2.36±0.07. The results of the best-fit models are reported in Tab. [5] and are shown in 
Fig. [T71 The spectral bias (see § |2j) is the strongest in the soft band while it is absent in the 
high-energy band, being already almost negligible above 1 GeV. 

From the logiV-logS' in the whole band we would expect (assuming a power law with a 
photon index of 2.4 and that the blazar population is not changing dramatically with energy) 
to find breaks at: 6.7xl0 -8 , 2.6xl0 -9 , and lxlO -10 ph cm -2 s _1 for the soft, medium, and 
hard bands respectively. Indeed these expectations are confirmed by the ML fits in the 
individual bands (e.g. see Tab. The hard band constitutes the only exception where the 
flux distribution barely extends below the flux at which the break might be observed. 

The average spectral properties of the sample change with energy. We find that the in- 
trinsic index, distribution is compatible with a Gaussian distribution with means of 2.25±0.02, 
2.43±0.03, and 2.17±0.05. In the three bands the fraction of BL Lac-to-FSRQ is: 0.61, 1.14, 
and 3.53 respectively with identification incompletenesses of 0.18, 0.25, and 0.25 respectively. 
It is apparent that the hardest band is the best one for studying BL Lac objects since the 
contamination due to FSRQs is rather small. 
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Fig. 17. — Source count distributions for the soft (0.1-1.0 GeV, left), medium (1.0-10.0 GeV, 
center) and high energy (10.0-100.0 GeV, right) band reconstructed with the method re- 
ported in § 16.21 

8. Contribution of Sources to the Diffuse Background 

The source count distribution can be used to estimate the contribution of point-like 
sources to the EGB emission. This allows us to determine the fraction of the GeV diffuse 
emission that arises from point-like source populations measured by Fermi. As specified in 
§ El this estimate does not include the contribution of sources which have been directly de- 
tected by Fermi since these are not considered in the measurement of the diffuse background. 
This estimate includes all those sources which, because the detection efficiency changes with 
flux, photon index and position in the sky, have not been detected. 

The diffuse emission arising from a class of sources can be determined as: 

where fi max is the geometrical sky area and the (1 — fl(T, S)/Q max ) term takes into 
account that the threshold at which LAT detects sources depends on both the photon index 
and the source flux. We note that neglecting the dependence of Q on the photon index (i.e. 
using the mono- dimensional sky coverage reported in Fig. [7]) would result in an underestimate 
of the diffuse flux resolved by Fermi into point-sources. The limits of integration of Eq. [12] 
are r min = 1.0, r max = 3.5, and S mSuX = 10~ 3 ph cm -2 s -1 . We also note that the integrand 
of Eq. [12] goes to zero for bright fluxes or for photon indices which are either very small 
or very large; thus the integration is almost independent of the parameters reported above. 
The integration is not independent of the value of S m i n which is set to the flux of the faintest 
source detected in the sample. For the analysis of the whole band S' m i n =9.36xl0 _10 ph 
cm -2 s _1 while for the low, medium and hard band S m j n is set to: 5.17xl0~ 9 ph cm -2 s _1 , 
3.58xl0~ 10 ph cm -2 s _1 , and 6.11xl0~ n ph cm -2 s -1 respectively. 
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Since in the measurement of lAbdo et al.l (j2010dl ) the sources which are subtracted are 
those detected in 9 months of operation, the coverage used in Eq. [12] is the one corresponding 
to the 9 months survey. The uncertainties on the diffuse flux have been computed by per- 
forming a bootstrap analysis. Integrating Eq. [TJ]we find that the point source contribution is 
1.63(±0.18) x 10 _6 ph cm~ 2 s _1 sr -1 where the systematic uncertainty is 0.6xl0~ 6 ph cm -2 
s -i sr -i This corresponds to 16( ±1.8) % (±7 % syst ematic uncertainty) of the Isotropic 
diffuse emission measured by LAT ( lAbdo et al.ll2010dl ) above 100 MeV. This small fraction 
is a natural consequence of the break of the source counts distribution. However, it is also 
possible to show that the parameter space for the faint-end slope 02 is rather limited and 
that a break must exist in the range of fluxes spanned by this analysis. Indeed, for a given 
02 (and all the other parameters of the log Af-log S fixed at their best-fit values) one can 
solve Eq. [12] to determine the flux at which the integrated emission of point sources exceeds 
the one of the EGB. Repeating this exercise for many different values of the 02 parameter 
yields an exclusion region which constrains the behavior of the log AMog S at low fluxes. 
The results of this exercise are shown in Fig. [16] From this Figure it is apparent that the log 
A"-log S must break between Fioo ~ 2 x 10 _9 ph cm 2 s _1 and Fioo ~ 6.6 x 10 _8 ph cm 2 s _1 . 
For a small break (e.g. 0\ — 02 ~ 0.2 — 0.3 and then 2 ~2.2-2.3), the integrated emission 
of point sources would still match the intensity of the diffuse background at Fioo ~ 10~ 9 ph 
cm 2 s -1 which are sampled by Fermi. Thus not only the break has to exist, but this simple 
analysis shows that it has to be strong (see also § 18 .3[) . not to exceed the intensity of the 
diffuse emission. 



The logAMogS 1 in the whole band goes deeper than the source count distributions 
derived in the smaller bands. This is clearly shown in Fig. [18] Given the fact that most 
of the source flux is emitted below 1 GeV (for reasonable photon indices) , the source count 
distribution in the soft band (0.1-1.0 GeV) is the one which gets closer to the logAMogS' in 
the whole band in terms of resolved diffuse flux. 



The log AMog S in the whole bands shows a strong break with a faint-end slope (e.g. 02) 
robustly constrained to be <2. In this case the integral reported in Eq. [T2l converges for small 
fluxes and it can be evaluated at zero flux to assess the maximum contribution of Fermi- 
like sources to the diffuse background. This turns out to be 2.39(±0.48) x 10 _6 ph cm -2 



s 1 sr 1 (1.26x10 6 ph cm 2 s 1 sr 1 systematic uncertainty) which represents 2 3(±5) 



(12% systematic uncertainty) of the Fermi diffuse background (lAbdo et al.ll2010dl ). This is 
a correct result as long as the log AMog S of point-sources (i.e. blazars) does not become 
steeper at fluxes below the ones currently sampled by Fermi. A given source population 
normally exhibits a source count distribution wi th a single downwards break (e.g. see the 
case of radio-quiet AGN in ICappelluti et al.ll2007l ). This break is of cosmological origin since 
it coincides with the change of sign of the evolution of that population. As can be clearly seen 
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in the redshift distribution in lAbdo et al.l ( l2010eh the epoch of maximum growth of blazars 



corresponds t o redshift 1.5-2.0 which co incides well with the peak of the star formation in the 



Universe (e.g. iHopkins &: Beacomll2006[ ). Since Fermi is already sampling this population it 
is reasonable to expect no other breaks in the source count distribution of blazars. Under this 
assumption, the result of the integration of Eq. [12] are correct. The results of this exercise are 
shown in Fig. [19] and summarized in Tab. [6] Since the 10-100 GeV source counts distribution 
does not show a break, its integral diverges for small fluxes. Thus, in both Fig. [19] and Tab. [6] 
we decided to adopt, as a lower limit to the contribution of sources to the diffuse emission 
in this band, the value of the integral evaluated at the flux of the faintest detected source. 

The different levels of contribution to the diffuse background as a function of energy 
band might be the effect of the mixing of the two blazar populations. In other words, as 
shown in § [7] FSRQs are the dominant population below 1 GeV while BL Lacs are the 
dominant one above 10 GeV. Given also that FSRQs are softer than BL Lacs (see also § [3]), 
it is naturally to expect a modulation in the blazar contribution to the diffuse emission as 
a function of energy. This can clearly be seen in Fig. [20] which shows the contribution of 
FSRQs and BL Lacs to the diffuse emission. This has been computed integrating the source 
count distribution of Tab. H] to the minimum detected source flux which is 9.36xl0 _10 ph 
cm -2 s" 1 and and l.llxl0 _8 ph cm -2 s _1 for BL Lacs and FSRQs respectively. It is clear 
that FSRQs are contributing most of the blazar diffuse emission below 1 GeV while BL Lacs, 
given their hard spectra, dominate above a few GeVs. The spectrum of the diffuse emission 
arising from the blazar class is thus curved, being soft at low energy (e.g. below 1 GeV) and 
hard at high energy (above 10 GeV), in agreement with the results of the analysis of the 
source count distributions in different bands. 



8.1. Additional Tests 
8.2. Source Count Distribution above 300 MeV 

The effective area of the LAT decreases quickly below 300 MeV while at the same time 



both the PSF size and the intensity of the diffuse background increase (e.g. see lAtwood et al. 



20091 ). In particular at the lowest energies, systematic uncertainties in the instrument re- 
sponse might compromise the result of the maximum likelihood fit to a given source (or set 
of sources). In order to overcome this limitation we constructed, with the method outlined 
in § El the log N~log S of point sources in the 300 MeV-100 GeV band. Considering that 
in the E> 100 MeV band the log N-log S shows a break around 6-7xl0~ 8 ph cm -2 s _1 and 
assuming a power law with a photon index of 2.4, we would expect to detect a break in the 
(E>300MeV) log AMog S around ~1.5xl0~ 8 ph cm -2 s -1 . Indeed, as shown in Fig. EU the 
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Fig. 18. — Contribution of point-sources to the diffuse GeV background. The red solid line 
was derived from the study of the logAMogS* in the whole band while the blue solid lines 
come from the study of individual energy bands (see §[7|). The bands (grey solid and hatched 
blue) show the total (statistical plus systematic) uncertainty. 



break is detected at 1.68(±0.33) x 10~ 8 ph cm" 2 s" 1 . Moreover, as Fig. |2"T1 shows, the break 
of the log iV-log S and the one of the sky coverage are at different fluxes. More precisely, 
the source counts start to bend down before the sky coverage does it. This is an additional 
confirmation, along with the results of § that the break of the log iV-log S is not caused by 
the sky coverage. The parameters of this additional source count distribution are reported 
for reference in Tab. 



8.3. Simulating a log N— log S without a break 



In order to rule out the hypothesis that the sources detected by Fermi produce most 
of the GeV diffuse emission, we performed an additional simulation. In this exercise the 
input log iV-log S is compatible with a single power law with a differen tial slope of 2 .23. At 
bright fluxes this log iV-log S is compatible with the one reported in lAbdo et al.l (j2009al ) 
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Fig. 19. — Contribution of point-sources to the diffuse GeV background obtained by extrap- 
olating and integrating the log iV-log S to zero flux. The red solid line was derived from the 
study of the logA^-logS* in the whole band while the blue solid lines come from the study 
of individual energy bands (see § [7]). The bands (grey solid and hatched blue) show the 
total (statistical plus systematic) uncertainty. The arrow indicates the lower limit on the 
integration of Eq. [12] for the 10-100 GeV band. 



and at fluxes F 10 o > 10 ph cm s accounts for ~70 % of the EGB. In this scenario the 
surface density of sources at F 10 o > 10~ 9 ph cm -2 s _1 is 0.8 deg -2 (while the one we derived 
in § 16.81 is 0.12 deg -2 ). To this simulation we applied the same analysis steps used for both 
the real data and the simulations analyzed in § HI Fig. [22] compares the flux distribution of 
the sources detected in this simulation with the distribution of the real sources detected by 
LAT and also with the sources detected in one of the simulations used in § [H It is apparent 
that the flux distribution of the sources detected in the simulation under study here is very 
different from the other two. 

Indeed, in the case point-like sources produce most of the EGB Fermi should detect 
many more medium-bright sources than are actually seen. A Kolmogorv-Smirnov test yields 
that the probability that the flux distribution (arising from the log AMogS' tested in this 
section) comes from the same parent population as the real data is < 10 -5 . This probability 
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Fig. 20. — Contributions of different classes of blazars to the diffuse GeV background ob- 
tained by integrating the log iV-log S. The red and the blues solid lines show the contribution 
of FSRQs and BL Lacs respectively, while the pink solid line shows the sum of the two. The 
bands around each line show the total (statistical plus systematic) uncertainty. 

becomes 5 x 10~ 4 if the x 2 t es t is used. The KS test between the flux distribution of one of 
the simulations used in § H] and the real data yields a probability of ~87 % that both come 
from the same parent population while it is ~91 % if the \ 2 test is used. 

Thus the hypothesis that Fermi is resolving (for F 10 o > 10~ 9 ph cm -2 s" 1 ) the majority 
of the diffuse background can be ruled out at high confidence. 

9. Discussion and Conclusions 

Fermi provides a huge leap in sensitivity for the study of the 7-ray sky with respect 
its predecessor EGRET. This work focuses on the global intrinsic properties of the source 
population detected by Fermi at high Galactic latitudes. 

We constructed the source count distribution of all sources detected above |6| >20°. 
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Fig. 21. — Source count distribution of all (TS> 25, |6| > 10°) sources in the 300 Me V- 
100 GeV band. The distribution has been multiplied by (F 100 /10~ 8 ) L5 . The dashed line 
shows the sky coverage (scaled by an arbitrary factor) used to derive the source counts. 
Note that the break of the log N - log S and that one of the sky coverage are at different 
fluxes. 

This distribution extends over three decades in flux and is compatible at bright fluxes (e.g. 
Fioo > 6 x 10~ 8 ph cm -2 s _1 ) with a Euclidean function. Several methods have been employed 
to show that at fainter fluxes the log iV-log S displays a significant flattening. We believe that 
this flattening has a cosmological origin and is due to the fact that Fermi is already sampling, 
with good accuracy, the part of the luminosity function which shows negative evolution (i.e. 
a decrease of the space density of sources with increasing redshift). This is the first time 
that such flattening has been found in the source count distributions of 7-ray sources and 
blazars. We also showed that the log iV-log S of blazars follows closely that of point source, 
showing that most of the unassociated high-latitude sources in the 1FLG catalog are likely 
to be blazars. At the fluxes currently sampled by Fermi (e.g. Fioo > 10 _9 ph cm -2 s _1 ) the 
surface density of blazars is 0.12^q;q2 deg -2 and this is found to be a factor ~4 larger than 
previous estimates. 

The average intrinsic spectrum of blazars is in remarkably good agreement with the 
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Fig. 22. — Flux distributions of detected sources (TS>50 and \b\ >20°) for three different 
realizations of the 7-ray sky. The solid thick line corresponds to a log N - log S distribution 
which resolves approximately ~70% of the GeV diffuse background, while the dashed line 
corresponds to the log iV-log S derived in this work which resolves approximately ~23 % of 
the diffuse background. For comparison the thin solid line shows the flux distributions of 
the real sample of sources detected by Fermi. 



spectrum of the GeV diffuse emission recently measured by Fermi (lAbdo et al.ll2010dl ). Nev- 
ertheless, integrating the log iV-log S, to the minimum detected source flux, shows that at 
least 16.0^!g % (the systematic uncertainty is an additional 7 %) of the GeV background can 
be accounted for by source populations measured by Fermi. This is a small fraction of the 
total intensity and it is bound not to increase dramatically unless the log N-logS becomes 
steeper at fluxes below 10~ 9 ph cm~ 2 s _1 . This generally never happens unless a different 
source class starts to be detected in large numbers at fainter fluxes. 



Thompson et al.l (120071 ) predict the integrated emission of starburst galaxies to be 
10~ 6 ph cm -2 s _1 sr _1 (above 100 MeV). This would represent ~10% of the LAT diffuse 
background and would be comparable (although a bit less) to that of blazars found here. 
Indeed, their prediction that M82 and NGC 253 wou ld be the first two starburst galaxies 
to be detected has been fulfilled (lAbdo et al.l l2010ai ). A similar contribution to the GeV 
diffuse background shoul d arise from the integrated emission of normal star forming galaxies 
(IPavlidou fc Fieldsll2002l ). In both cases (normal and starburst galaxi es) 7-ravs are produ ced 
from the interaction of cosmic rays with the interstellar gas (e.g. see lAbdo et al.ll2009bl ). It 
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is natural to expect that both normal and starburst galaxies produc e a fraction of the d iffuse 
emission since now both classes are certified 7-ray sources (see e.g. lAbdo et al.ll2010bl ). 



It is also interesting to note that pulsars rep resent the second largest popul ation in 
our high-latitude sample (see Tab. [2]). According to iFaucher-Giguere fc Loebl (120091 ) pulsars 
and in particular millisecond pulsars can produce a relevant fraction of the GeV diffuse 
emission. However , given the strong break, typically at a few GeVs, in their spectra (e.g. see 
Abdo et al.ll2010ff ). millisecond pulsars are not expected to contribute much of the diffuse 
emission above a few GeVs. Finally radio-quiet AGN might also contribute to the GeV 
diffuse background. In these objects the 7-ray emission is supposedly pr oduced by a non- 



thermal elect rons present in the cor ona above the accretion disk (see e.g. Ilnoue et al.l 12008 



for details). Ilnoue &: Tot anil ( 120091 ) predict that, at fluxes of F 100 < 10 °ph cm 



radio-quiet AGN outnumber the blazars. According to their prediction, most of background 
could be explained in terms of AGN (radio-quiet and radio-loud). 

It is thus clear that standard astrophysical scenarios can be invoked to explain the GeV 
extragalactic diffuse background. However, the main result of this analysis is that blazars 
account only for <40 % of it0- It remains a mystery why the average spectrum of blazars is so 
similar to the EGB spectrum. Taken by itself, this finding would lead to believe that blazars 
might account for the entire GeV diffuse background. However, we showed (see Fig. [22] and 
§ 18.31 for details) that in this case Fermi should have detected a much larger number (up 
to ~50%) of medium-bright sources with a typical flux of F 100 > 10~ 8 ph cm -2 s" 1 . This 
scenario can thus be excluded with confidence. Thus, the integrated emission from other 
source classes should still have a spectrum declining as a power-law with an index of ~ 2.4. 
This does not seem to be a difficult problem to overco me. Indeed, at least in t he case of star 



formin g galaxies we note that in the modeling of both iFields et al.l (120101 ) and iMakiya et al. 



( 120101 ) the integrated emission from these sources displays a spectrum similar to the EGB 
one (at least for energies above 200 MeV). Moreover, in this work we also found that the 
contribution to the diffuse emission of FSRQs and BL Lacs is different, FSRQs being softer 
than BL Lacs. Thus, the summed spectrum of their integrated diffuse emission is curved, 
softer at low energy and hard at high (> 10 GeV) energy. This makes it slightly different from 
the featureless power-law of the diffuse background. All the estimates presented here will 
be refined with the derivation of the blazar luminosity function which is left to a follow-up 
paper. 



4 This includes extrapolating the source counts distribution to zero flux and taking into account statistical 
and systematic uncertainties. 
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Table 3. Results of the power-law fits to the differential source count distributions 
obtained with the standard method of § 16. f I 







Sample Limits 




Best-fit Parameters 




SAMPLE 


# Objects 


TS> 


|b|> 


A a 


ft 


s b b 




ALL 
ALL 


425 
483 


50 
50 


20° 
15° 


i ic+0.15 
i - io -0.15 

i 74+016 


2 63+ 22 
2.60±^ 9 7 


6-40±i;°| 


1 fi4+0-06 

1 fi0+ 006 



a In units of 10 14 cm 2 s deg 2 . 

b In units of 10~ 8 ph cm~ 2 s" 1 (0.1<E<100 GeV). 



Table 4. Results of the best fits to the source count distributions. 



Sample Limits Best-fit Parameters 



SAMPLE 


# Objects 


Incompl. 


TS> 


|b|> 


A a 


Pi 


s 6 b 


ft 


M 


a 


ALL 


425 





50 


20° 


16.46±0.80 


2.49±0.12 


6.60±0.91 


1.58±0.08 


2.36±0.02 


0.27±0.01 


BLAZAR 


352 


0.13 


50 


20° 


18.28±1.00 


2.48±0.13 


7.39±1.01 


1.57±0.09 


2.37±0.02 


0.28±0.01 


FSRQ 


161 


0.19 


50 


20° 


72.41±5.76 


2.41±0.16 


6.12±1.30 


0.70±0.30 


2.48±0.02 


0.18±0.01 


BL Lac 


163 


0.19 


50 


20° 


0.106±0.009 


2.74±0.30 


6.77±1.30 


1.72±0.14 


2.18±0.02 


0.23±0.01 


Unassociated 


56 





50 


20° 


3.12(±0.5) x 10~ 5 


3.16±0.50 


4.48±1.3 


1.63±0.24 


2.29±0.03 


0.20 



a In units of 10 14 cm 2 s deg 2 . 
b In units of 10 -8 ph cm -2 s _1 . 



Table 5. Results of the best fits to the source counts distributions in different Energy bands. Parameters without 

error estimate were kept fixed during the fitting stage. 









Sample 


Limits 






Best-fit Parameters 






BAND 


# Objects 


Incompl. 


TS> 


|b|> 


A a 


fh. 




02 


/* 




0.1-1.0 GeV 
1.0-10.0 GeV 
10.0-100.0 GeV 
0.3-100.0 GeV 


362 
597 
200 
759 








25 
25 
25 
25 




o o o o 


4.00±0.21 
1.097±0.05 
8.3(±0.6)xl0- 3 
5.33±0.19 


2 5^+ ' 17 
Z.OO_ 22 

2 44+ - 15 


- 7C :+0.44 
,J -'' J -2.22 

0.23±0.06 

1 69+ ' 33 


, oo+0.13 
J--3»_ . 46 

1.52±?;f 

i 70 +0.06 
' u -0.07 


9 9C-+0.02 
z - ZJ -0.()2 

2.43 
2.17±0.05 

9 or+0.02 


32+ 01 

0.40tg;g 

0.82±°;°° 

o.sotHi 



a In units of 10 14 cm 2 s deg 2 . 
b In units of 10~ 8 ph cm~ 2 s _1 . 
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Table 6. Diffuse emission arising from point sources. The lower part of the table shows 
the values of the integrated emission when the source counts distributions are extrapolated 
to zero flux. Errors are statiscal only (see § [5] for a discussion about systematic 

uncertainties). 



Band 

(GeV) 


EGB Intensity 11 

(ph cm -2 s — 1 sr _1 ) 


Point Source Diffuse Emission 

(ph cm -2 s — 1 sr _1 ) 


Fraction of EGB Intensity 

(%) 


a . b 

(ph cm -2 s^" 1 ) 


0.1-100 
0.1-1.0 
1.0-10 
10-100 


1.03xl0" 5 
9.89xl0~ 6 
3.85xl0~ 7 
1.50xl0~ 8 


1.63(±0.18) x 10~ 6 
1.54±°- 2 J x 10- 6 
2.93+J-?? x io- 8 
1.36±g;|t x 10- 9 


15.8(±1.8) 
15.5t?- 9 3 
7 6 +5 -° 

' - D -1.8 

9.0l|J 


9.36 
51.1 
3.58 
0.61 


0.1-100 
0.1-1.0 
1.0-10 
10-100 


1.03xl0~ 5 
9.89xl0~ 6 
3.85xl0~ 7 
1.50xl0~ 8 


2.39(±0.48) x 10~ 6 
2.07±g;|f x 10~ 6 
5-49±f;?g x IO" 8 
> 1.36 x 10~ 6 


22.5(±1.8) 
20.9™ 
14.2t^ 4 2 
> 9.0 


o o o o 



a The intentisities of the EGB emission are derived from lAbdo et al.l (|2010dh . 

b Lower flux of integration of the source counts distributions (see Eq. [T2J in units of 10 _10 ph cm -2 s _1 . 

c The source counts distribution in the 10-100 GeV does not show a break and thus, its integral to zero flux 
diverges. As a lower limit on the diffuse emission, we adopted the value computed at the faintest detected 
source. 
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